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Given the promising features of the recently proposed Barcelona-Catania-Paris (BCP) functional 
1], it is the purpose of this paper to still improve on it. It is, for instance, shown that the number of 
open parameters can be reduced from 4-5 to 2-3, i.e. by practically a factor of two. One parameter 
is tightly fixed by a fine-tuning of the bulk, a second by the surface energy. The third is the strength 
of the spin-orbit potential on which the final result does not depend within the scatter of the values 
used in Skyrme and Gogny like functionals. An energy rms value of 1.58 MeV is obtained from a 
fit of these three parameters to the 579 measured masses reported in the Audi and Waspra 2003 
compilation. This rms value compares favorably with the one obtained using other successful mean 
field theories. Charge radii are also well reproduced when compared with experiment. The energies 
of some excited states, mostly the isoscalar giant monopole resonances, are studied within this model 
as well. 



I. INTRODUCTION 

The physics of atomic nuclei and nuclear systems is 
very rich but extremely complicated to describe. It 
unites many of the toughest aspects of the quantum many 
body problem (common to other branches of physics, 
like condensed matter physics) in one spot. To name 
only a few, we can mention the existence of four dif- 
ferent types of fermions (proton/neutron-spin up/down), 
the relevance of pairing and quartet correlations in the 
dynamics, or the role of the mean field approach being 
only a first approximation while quantum fluctuations 
are very important as a consequence of the smallness of 
nuclei. Moreover, the nucleon-nucleon (NN) force is very 
complicated with its hard core, tensor (dipolar) compo- 
nents and strong spin orbit interactions. It is remarkable 
that nuclear theory has nowadays mastered a great deal 
of those challenges and many of the experimental data 
are very well described either on a phenomenological ba- 
sis, or, more rarely, with microscopic approaches. In this 
work we try to progress with the latter strategy. 

Theories aiming towards a global description of low en- 
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ergy nuclear data always start with a mean field descrip- 
tion and a semi-phenomenological nuclear energy density 
functional. In a recent paper we proposed an ap- 
proach slightly different from the usual for the establish- 
ment of a nuclear energy density functional. We tried to 
follow as much as possible the Kohn Sham Density Func- 
tional Theory (KS-DFT) which is in common use, e.g., in 
atomic and molecular physic 043 ■ The basis of this the- 
ory lies in the Hohenberg-Kolm theorem [8j , which states 
that for a Fermi system, with a non-degenerate ground 
state, the total energy can be expressed as a functional 
of the local density p(r) only. Such a functional reaches 
its variational minimun when evaluated with the exact 
ground state density. Unfortunately, this is an existence 
theorem and does not include a constructive algorithm 
to build the energy density functional. The usual ap- 
proach is provided by the KS method @ which intro- 
duces a set of auxiliary single-particle states and takes 
for the kinetic energy the Slater form in this basis. In 
the KS method the energy density functional is written 
as a sum of two terms. One of them corresponds to the 
uncorrelated kinetic energy density plus the direct contri- 
bution derived from the underlying two body (Coulomb) 
interaction. The other piece contains an unknown term 
for the exchange and correlation energies. For a practi- 
cal implementation of the method, the latter term has 
to be guessed and therefore many strategies have been 
developed along the years to describe different physical 
phenomena. A popular strategy in atomic and molec- 
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ular physics is to use accurate theoretical calculations 
(mainly by Monte-Carlo techniques) as a function of the 
density in simplified models where finite size effects are 
neglected. The results of those calculations are used to 
make an educated guess of the unknown exchange and 
correlation part of the functional in the real system. In 
this way the minimization procedure gives Hartree-like 
equations for the single-particle states , where the poten- 
tial part includes in an effective way the overall exchange 
and correlation contribution. 

The application of the KS-DFT scheme to self bound 
systems like nuclei, is not straightforward. This hinges 
on the fact that, contrary to atomic or molecular systems 
where the ions form a natural center, in nuclei mean field 
type of equations always break some symmetries as, for 
instance, translational and rotational invariance, but also 
eventually particle number, etc. This problem has, how- 
ever, recently been solved very elegantly in Ref. |10l - [l2j 
where it is shown that the Hohenberg-Kohn theorem [8| , 
on which the KS-DFT scheme is based, has to be refor- 
mulated for the 'intrinsic' system in the case of nuclei or 
other self bound quantum droplets. 

In nuclear physics to build up the exchange and corre- 
lation contribution amounts to a detailed determination 
of the ground state of symmetric nuclear matter and of 
neutron matter to be used as a guide for finite nuclei en- 
ergy density functionals. Accurate nuclear matter calcu- 
lations are very complicated and not as much advanced as 
for electron systems of condensed matter. Nonetheless, 
in recent years, quite some progress has been achieved 
and we will build our KS-DFT approach on the equation 
of state (EOS) of Baldo et al., developed in Refs. (l3l — 
fiol ] , which is based on the well known hole line expansion. 
This EOS is calculated from realistic two body and three 
body interactions treated at the Brueckner Hartree-Fock 
level. It is compatible with phenomenological constraints 
coming from data of heavy ion reactions [l9j as well as 
from the analysis of astrophysical observations [2(j. It 
is worth mentioning that the only work which tried to 
follow the same strategy as here, is the one by Fayans 
[2l| . However, that work was based on a somewhat older 
EOS [22] and had remained an isolated work. 

In the present work, we will follow along the same lines 
as in pj, i.e. with a Kohn-Sham approach based on the 
same microscopic bulk input. However, we will be able 
to reduce the number of parameters by two without los- 
ing accuracy. As before, we will consider an additional 
bulk parameter for a precise adjustment of the numeri- 
cally obtained E/A value. However, for the surface, we 
will reduce the number of parameters from three to one 
by the condition that in infinite symmetric and infinite 
pure neutron matter the strength parameters of the finite 
range term reproduce the coefficients of the quadratic 
terms of the bulk energy density obtained from the mi- 
croscopic calculation. Therefore, no subtraction term in 
the finite range part as in [l[ is necessary. So, for the sur- 
face the range ro remains the only adjustable parameter, 
see the more detailed discussion in the main text. Then, 



we are left with only two adjustable parameters, one for 
the bulk and one for the surface. Although, in principle, 
there are, together with the strength of the spin-orbit 
term, three parameters, we want to emphasize that the 
main parameters are the ones from bulk and surface. We 
stress this point, because it seems to us a reduction to 
two basic physical inputs to the binding energy of nuclei: 
energy per particle of infinite matter and surface energy! 
That this is possible is as much surprising as it is gratify- 
ing! It also should be mentioned that the adjustment of 
those two parameters is extremely sensitive: both have 
to be fine tuned to the order of 10 -4 ! Any minor devi- 
ation from the nominal values very rapidly and strongly 
deteriorates the results! 

It should be pointed out, however, that KS-DFT ad- 
dresses only the ground state and is, in principle, not 
taylored to describe excitations of the system. Never- 
theless, we, tentatively, will also use it, at least, for the 
description of the Giant Monopole Resonances (OMR's) 
where the fact that in KS-DFT the effective mass is equal 
to the bare mass, i.e. m* = m, plays a minor role. 

The paper is organized as follows: The following sec- 
tion is devoted to briefly recall our previous BCP func- 
tional. In the third section, the new energy density func- 
tional built up in this paper is discussed. The results 
obtained with this improved functional, that we will call 
BCPM (Barcelona-Catania- Paris-Madrid), are also pre- 
sented in the same section. The predictive power of the 
BCPM functional regarding other observables such as 
charge radii, quadrupole and octupole deformations and 
fission barriers are discussed in the fourth section. The 
ability of our BCPM functional for describing for instance 
the isoscalar giant monopole resonance is the subject of 
the fifth section. Finally, the summary and conclusions 
are given in the last section. 



II. FORMER BCP FUNCTIONAL 

The former BCP functional was proposed in Ref. 
This and subsequent refinements [23] are based on the 
Kohn Sham Density Functional Theory (KS-DFT) where 
the one body density p(r) plays a central role. In the KS- 
DFT theory an auxiliary set of A orthonormal wave func- 
tions <fi(r), where A is the mass number, is introduced to 
express formally the density as if it were obtained from 
a Slater determinant as a sum of the product of single 
particle wave functions 

p(r)=5> 2 (r)| 2 , 

i 

with the ip's determined from the minimization of the 
ground state energy. In condensed matter and atomic 
physics the energy density functional E[p(r)] is usually 
split into two parts @: 

E = T [p] + W[p}. 
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The first piece To[p] corresponds to the uncorrelated ki- 
netic energy written in the usual manner as 

i 

The other piece, W[/?], contains the potential energy and 
the correlated part of the kinetic energy. In the BCP 
family of functionals this part is given as the sum of 
the spin-orbit energy, the Coulomb part and the nuclear 
energy term which depends on the neutron and proton 
densities p n and p p , respectively. In nuclear physics, con- 
trary to the situation in condensed matter and atomic 
physics, the contribution of the spin-orbit interaction to 
the energy functional is very important. Non-local con- 
tributions have been included in DFT in several ways 
long ago (see Ref. [25] for a review of this topic). In 
line with this, we split the spin-orbit part into an un- 
correlated part E s -°- plus a remainder. The form of the 
uncorrelated spin-orbit part is taken exactly as in the 
more modern Skyrme [H, H3| or Gogny forces [28L [HI • 
Another piece that we explicitly split off is the Coulomb 
part Ec- We treat this contribution at lowest order, i.e. 
the direct term plus the exchange contribution in the 
Slater approximation, that is 

E " = \ /^ r %(r)|r - r'|-Vp(r'), 

and 

^=-l(l) v 7^w v8 

with E c = Eg + E C X . 

The nuclear energy functional contribution E^ n i^p ni pp\ 
contains the nuclear potential energy as well as additional 
correlations. This contribution is divided into a finite 
range term Ef^\p ni p p \ to account for correct surface 
properties and a bulk correlation part E°% t [p n , p p ] that 
we take from one of the most advanced microsc opically 
determined EOS existing so far in the literature |13T - fl9j 
as mentioned before. Collecting all these contributions 
our final KS-DFT functional reads 

E = T + E s -°- + E? nt + E[J + E C . (1) 

One of the prominent results of our recent work within 
the KS-DFT BCP density functional scheme Q was that, 
with a very reduced number of adjustable parameters, 
rms values for binding energies and radii came out to be 
of the same quality as the ones of most of the more mod- 
ern Skyrme [H E3 or Gogny [H, [H| functionals. The 
reduction of the number of parameters stemmed from 
the fact that, prior to a fit of the functional to nuclear 
masses, the bulk part of the functional E°^ t was adjusted 
to the microscopic EOS obtained in Refs. [13l4l9| by fit- 
ting polynomials in the total density p = p n + p p to mi- 
croscopic results in symmetric and pure neutron matter, 



followed by a quadratic interpolation in the asymmetry 
parameter j3 — (p n — p p )/p between these two limits. It 
has to be stressed that the use of a polynomial in density 
obeys to practical reason. This fit procedure necessarily 
has a slight freedom, since the numerical grid points are 
not dense and not 100 % fixed, neither from the theo- 
retical nor from the numerical side. We will exploit this 
freedom, which is of the order of 10 of the energy per 
particle E/A, to fine tune our results. Once the bulk 
part is adjusted in this way, it will not be varied any 
more in the parameter search which determines the nu- 
clear masses by adding to the functional a finite range 
term. 

The additional surface term in [1] was of the form 

Ef n ?[Pn,P P ] = \ W/drdrV t (r)« M ,(r-r> t >'X2) 

M' 

1 t, t > J 

where the index t is the label for neutron and proton, 
i.e. t = n,p, and jt,t' the volume integral of Vt.t'ij"). The 
subtraction term was introduced not to contaminate the 
bulk part, determined from the microscopic infinite mat- 
ter calculation. For the finite range form factor Vtfi ( r ) a 
simple Gaussian ansatz 

v t ,Ar) = V t ^e- r2 / r « 

was made. The strength parameters were chosen to dis- 
tinguish only between like and unlike particles, i.e. 

Vp,p = V n , n = Vl; V n , p = V P: n = Vu, 

so that the finite range term contained three adjustable 
parameters: Vl,Vu, and tq- Together with the adjust- 
ment of the bulk and the spin-orbit strength, the func- 
tional in [l[ has all in all five adjustable parameters, even 
though they are not all to be considered on the same level 
of significance. We will come back to these considerations 
later in the text where we also will show how to avoid the 
subtraction term in the finite range part of the energy. 

For describing open shell nuclei one has to introduce 
pairing correlations. The formal way of including pair- 
ing within the DFT is also known since long from the 
generalization of the Hohenberg-Kohn theorem to paired 
systems In Q we have introduced pairing in our 

functional through the BCS approach using a zero-range 
density dependent interaction adjusted to reproduce the 
gap s of the Gogny force in neutron matter with m* = m 

We would like to point out that for other energy den- 
sity functionals, like the ones of the Skyrme, Gogny, or 
relativistic mean field type, the number of parameters 
seems much higher, typically more than ten. However, 
many of those parameters are implicitly used to get a 
reasonable nuclear matter and neutron matter equation 
of state. The advantage of the KS-DFT procedure is 
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that one clearly separates the tasks of reproducing the 
nuclear matter EOS from the most prominent finite size 
effects, namely the surface and Coulomb energies and the 
spin-orbit potential. 

The free parameters of the energy density functional 
reported in [l| were fitted either to reproduce the bind- 
ing energies of a few spherical nuclei or binding ener- 
gies and charge radii of the same set of spherical nuclei. 
We called these functionals BCP1 or BCP2 (Barcelona- 
Catania-Paris), respectively. The predictive power of 
these functionals was illustrated by computing the bind- 
ing energies and charge radii of 161 spherical nuclei be- 
tween 16 Ne and 224 U. We found that the rms in energy 
and charge radii are similar to the ones obtained with suc- 
cessful mean-field models such us the Gogny D1S force 
[2^ | , the Skyrme SLy4 interaction [2(1 and the relativistic 
mean field parametrization NL3 [32]. 

In subsequent publications we explored the proper- 
ties of BCP1 and BCP2 in describing quadrupole and 
octupole deformation ground state properties, fission, 
excited octupole states, etc (33l - [35j . Also some other 
parametrization of the functional based on fitting to the 
binding energy of deformed nuclei instead of the fitting 
to spherical nuclei, have been considered |23j. 



III. IMPROVEMENTS OVER BCP AND 
RESULTS 

A. New polynomial fitting 

Following [l|, [23| we write the bulk part of the energy 
density functional as 

E£t\f>P,Pn] = J dr[P s (p)(l - P 2 ) + P n (p)f3 2 ] P , (3) 

where the interpolating polynomials for symmetric and 
pure neutron matter, P s (p) and P n (p) respectively, read 



p s(p) = y2 a n[ — 



P n (p) = J2 b n( — ) > ( 4 ) 



where the reference densities are po=0.16 fm -3 and 
Pon— 0.155 fm -3 , respectively. 

As compared to the former version of the BCP energy 
density functional, the fit of the microscopic EOS has 
been redone in order to avoid some rather strong oscil- 
latory behavior in the density-dependent incompressibil- 
ity of the nuclear matter (see Section IV below) when 
plotted as a function of density. To cure this unwanted 
effect, a unique five order polynomial in the total den- 
sity is chosen for fitting the microscopic EOS in sym- 
metric nuclear matter in a wider range of densities, up 
to 0.625 fm -3 , instead of the prescription given in [l|, 
see also (23j. Furthermore special care has been paid to 
the smoothness not only of the fitted EOS but also of 
the density-dependent incompressibility. This has been 



n 


bn 


On(E/A=16) 


a n (E/A = 15.98) 


1 


-34.9726153509151 


-73.2920258 


-73.3826730960756 


2 


22.1823070966258 


49.9649055 


50.2977981668080 


3 


-7.15175564221894 


-18.0376009 


-18.3667336471822 


4 


1.79087375287097 


3.48617639 


3.60835940133109 


5 


-0.169591425337085 


-0.243552099 


-0.258847229282364 



Table I: Parameters (in MeV) of the polynomial fits, P s (p) = 
^2 n Om(p/po) n and P n {p) = J2 n b n (p/pon) n for symmetric 
and neutron matter, respectively. The reference densities are 
po = 0.16 fm -3 and po n = 0.155 fm -3 . Two sets of parame- 
ters for symmetric matter are given, they correspond to EOS 
with a minimum at the given values of E/A. 



achieved by keeping under control high order derivatives 
of the physical quantities. 

Since we want to construct the KS-DFT functional, as 
much as possible, on the basis of the microscopic calcula- 
tions, the bulk part E°^ t of the functional, directly related 
to bare NN and TBF (three-body forces), is determined 
once and for all and we will use it in ([3]) together with the 
local density approximation. However, as we mentioned 
in the Introduction, fine tuning of the fitting with differ- 
ent saturation points are performed in order to optimize 
the performance of the functional. We explored the in- 
terval between 15.97 and 16.03 MeV for the saturation 
energy, having fixed the saturation density at 0.16 fm -3 . 
The optimal value for the saturation energy per particle 
turns out to be 15.98 MeV. 

The values of the coefficients of the fitted polynomi- 
als for neutron matter and for two choices of symmetric 
matter leading to slightly different values of E/A are re- 
ported in table U A nice feature of the coefficients a n is 
that they follow a perfect linear dependence when plot- 
ted as a function of E/A and therefore the two sets of a n 
values are enough to reproduce the whole range of E/A 
values considered with only one parameter. 

As an aside, it may be important to mention that the 
density dependent part of the functional depends on in- 
teger powers of the density only. This is advantageous 
when trying to overcome the self-energy problem that 
plagues the use of theories beyond the mean field based 
on standard nuclear energy density functionals (36l |37| . 

As discussed in [38[ , the low density behavior of the mi- 
croscopic nuclear matter EOS has a characteristic trend, 
usually not reproduced by Skyrme and Gogny functionals 
( see also ref. |18l. l23j). missing there quite a substantial 
part of binding. We show the microscopic EOS for nu- 
clear and neutron matter as well as the corresponding 
polynomial fits in Fig. [T] In the lower panel of the same 
figure, we also display the symmetry energy computed 
using the polynomial fits. 
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later use, it will also be convenient to introduce the coef- 
ficient K' connected to the so-called skewness coefficient 
H3 by 



K' 



-Q = -27 p 



d\H/ P ) 
dp 3 i 



(7) 



Other relevant quantities in asymmetric nuclear matter 
are the symmetry energy and its first and second deriva- 
tives with respect to the density which govern the isovec- 
tor part of the nuclear interaction. The symmetry energy 
is defined as 



E 



1 d 2 fH 



sym 



(8) 




0.5 



density (fm ) 



Figure f : Upper panel: EOS of symmetric and neutron mat- 
ter obtained by the microscopic calculation (squares) and the 
corresponding polynomial fits (solid lines). For comparison 
the microscopic EOS of Refs. [24| are also displayed by open 
circles. Lower panel: Symmetry energy obtained from the 
polynomial fits 



B. Infinite matter properties 

The pressure P and density-dependent incompressibil- 
ity K in asymmetric nuclear matter are defined as (39| 



p = p2 d(H/ P ) m 



dp 



= n ' 

op 



and 



K = 9p 



d 2 H 
dp 2 



WP) 18 
9 dp 2 + p 



(5) 



(G) 



respectively, where H(p,/3) denotes the energy density 
in asymmetric nuclear matter. For symmetric mat- 
ter (/3 = 0) at saturation the K(p) defined by ((6]) re- 
duces to the well known incompressibility modulus Kq — 
9p 2 d \p=p ■ where po is the saturation density. For 



At saturation density one defines the symmetry energy 
coefficient as J — E sym (p ) and two coefficients more, L 
and K sym , directly related to the first and second deriva- 
tives of ([H]) with respect to the density at saturation, re- 
spectively 



L = 3 dEsyjn 

dp I 



2 ^' "sym 



dp 2 



\p=P 



(9) 



The values of the incompressibility modulus K$ and the 
coefficient K' defined in symmetric nuclear matter as well 
as the coefficients J, L and K sym which are related to the 
symmetry energy are displayed in Table HTl 

The binding energy and the saturation density in nu- 
clear matter are constrained by nuclear masses and elec- 
tron scattering experimental data. The range of ac- 
cepted values of the incompressibility modulus K are 
constrained by the experimental excitation energies of the 
isoscalar giant monopole resonance in finite nuclei since 
the pioneering works of Bohigas and collaborators (4l| 
and Blaizot [42|. However, different estimates of Kq us- 
ing different mean-field models predict slightly different 
values of this coefficient. Recently, a value Kq = 230 ±30 
MeV has been proposed [43] as a compromise among the 
different available estimates. The value of the symmetry 
energy coefficient J , that rules the isospin dependence 
of the nuclear interaction, is constrained by experimental 
data of heavy-ion collisions, pigmy dipole resonances and 
analog states. A range 30 < J < 35 MeV has been re- 
cently proposed for this coefficient [44| . The density con- 
tent of the symmetry energy is a relevant physical quan- 
tity related with many phenomena not only in terrestrial 
nuclei but also in neutron stars. Since the celebrated cor- 
relation established by Brown [45j between the slope of 
the symmetry energy and the neutron skin in 208 Pb, a 
considerable effort has been devoted to constrain the L 
parameter from available data. Antiprotonic atoms, nu- 
clear masses, heavy-ion collisions, giant and pigmy dipole 
resonances, proton -nucleus scattering as well as theo- 
retical calculations using a microscopic interactions have 
been used to extract the L coefficient (see Ref. [46[ and 
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references therein). From a compilation of the exist- 
ing data a range of values L = 55 ± 25 MeV has been 
proposed flfl Vu\ . The other two coefficients that char- 
acterize asymmetric nuclear matter around saturation, 
namely K' and K sym are more uncertain. An estimate 
K' = 700 ±500 MeV has been proposed in Ref. 0. The 
curvature of the symmetry energy, K sym , can be inferred 
from some non-relativistic [48( and relativistic [49j] calcu- 
lations. From these results a range —200 < K sym < 150 
MeV can be inferred. The nuclear matter properties of 
our BCPM energy density functional, that are obtained 
from a microscopic EOS calculation, lie within the ac- 
cepted ranges of values of the different quantities as it 
can be seen from a direct inspection of Table [Til 



Table II: Infinite nuclear matter properties of BCPM. All the 
parameters, except po (in fm ) and the dimensionless effec- 
tive mass m/m* , are given in MeV 



B/A po m/m* J 







K' K a 



-15.98 0.16 1.00 31.90 52.96 212.4 879.6 -96.75 



C. Surface term strengths derived from nuclear 
matter data and spin-orbit contribution 

In the new version of the functional the gaussian func- 
tions in the surface term can, in principle, have different 
ranges depending on the isospin channel tql and r$u in- 
troducing in this way a new parameter as compared to 
the old BCP1. However, the strengths Vl and Vjj are 
now no longer considered as free parameters. They have 
been determined in such a way that the bulk limit of the 
surface term in ([3]), that is 



/./' 



7m' / d d rp t (r)p tl ((r), 



reproduces the p 2 term of the bulk part of the energy 
density in Eq. ([3]). Therefore, no subtraction as in ^ 
is required. This procedure imposes the following re- 
lationships between the surface term strengths and the 
coefficients a\ and b\ of the polynomial fits of symmetric 
and neutron matter, respectively 



Vl 
V v 



2&i 



4ai 



r OLPO 

- 2bi 



(10) 
(11) 



Here tql and tqu are the ranges of the gaussians in the 
surface term, as defined above. The parameter b\ defined 
as bi = bi-^- has also been introduced to take into ac- 
count the different reference densities for symmetric and 
neutron matter. 
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Figure 2: Binding energy differences AB = B t h — £>Ex P (in 
MeV) as a function of neutron number JV for two choices of 
E/A, namely E/A — 15.97 in the upper panel and E/A = 
16.03 in the lower one 



On the other hand, the spin-orbit strength is practi- 
cally fixed from the outset because, within a certain mar- 
gin, the final result is nearly independent of it. We will 
take the spin-orbit strength from the value in the Gogny 
force with Wls = 130 MeV but, since we take the bare 
mass in the BCPM-functional, we have to reduce it by 
30 percent, i.e. we take in the present work the definite 
value of W LS = 90.5 MeV. 

From this discussion we see that the number of pa- 
rameters in the surface term gets reduced from three to 
two, namely the Tol and r^u ranges and that the spin- 
orbit strength Wls is practically fixed from the begin- 
ning. However, as we will see below, the actual fit allows 
for Tql = fou = r o and, thus, one is left with one param- 
eter only. 



D. Fitting protocol for finite nuclei 

Owing to the simplicity of the BCP energy density 
functional we are in the position to perform mass ta- 
ble evaluations on a personal computer in a reasonable 
amount of time ( a few hours ). Therefore we have de- 
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cided to include in the fit not only the masses of spher- 
ical nuclei but also the ones corresponding to deformed 
ground states as recently proposed The only limi- 

tation has been the restriction to even-even nuclei. Fi- 
nite nuclei properties are computed within the mean field 
framework provided by the well known Hartree-Fock- 
Bogoliubov (HFB) method [5l|. The pairing interaction 
is the same as in the original formulation of the functional 
[l[ but using the HFB instead of the BCS approach. 

Details of the HFB calculations are as follows: the 
quasiparticle operators have been expanded in a har- 
monic oscillator basis with the number of shells depend- 
ing on the nuclei considered. For nuclei with Z < 50 
we take eleven shells, for 50 < Z < 82 we consider 13 
shells and for Z >82 a 15 shell calculation is used. For 
the binging energy, an extrapolation to an infinite ba- 
sis is performed in the spirit of Ref. [52] . For each 
nucleus we compute the HFB energy with N and N+2 
shells and estimate the N — oo energy by the formula 
E(oo) = 2E(N + 2) - E(N). 

For the solution of the non-linear HFB equation, an 
approximate second order gradient method is used [53|. 
The information about the energy curvature reduces the 
number of iterations substantially as compared to other 
methods. 

As it is customary in this kind of calculations, the 
Coulomb exchange contribution is computed in the Slater 
approximation [54] and the Coulomb anti-pairing effect is 
not explicitly considered (see [55| for a discussion of this 
issue). The two body kinetic energy correction, which 
is typically considered as a way to correct for the lack 
of translational invariance of the whole procedure, has 
been taken into account with the pocket formula of Ref. 
[EH [57j as in previous versions of the BCP functional. 

In deformed nuclei the correlation energy stemming 
from symmetry restoration of the rotational invariance 
(the rotational energy correction) plays a relevan role. It 
rapidly changes from magic or semi-magic nuclei, where 
it is zero or very small, to strongly deformed mid-shell 
nuclei where it reaches values that can be as large as 6 
or 7 MeV in heavy nuclei. This is a consequence of its 
almost linear dependence with quadrupole deformation 
and therefore its inclusion is relevant to describe masses 
all over the Nuclide Chart. The correct way to com- 
pute this quantity is by evaluating the angular momen- 
tum projected energy associated with each intrinsic state. 
This is a tremendous task that can fortunately be alle- 
viated by considering an approximation to the projected 
energy that is obtained in the spirit of the Topological 
Gaussian Overlap Approximation. This formula, which 
is similar to the rotational energy correction, can be eas- 
ily evaluated at the mean field level (see Ref. (H for 
details). 

Finally, as already mentioned, for the pairing part we 
have taken the density dependent zero range force of Ref. 
[3l| that was devised to mimic the behavior of Gogny Dl 
pairing gap in nuclear matter. We have taken care of our 
effective mass equal to one by renormalizing the pairing 



force strength given in [3l| . 

We have used the 2003 evaluation of Audi and Wap- 
stra [59] for the experimental binding energies to be fitted 
and considered the set of 579 even-even nuclei where the 
binding energy has been measured (we exclude from our 
analysis the 193 extrapolated binding energies included 
in Audi and Wapstra's work, although we will use them 
to assess the quality of our fit). It has to be emphasized 
that nothing else but the rms deviation for the binding 
energy a(E) has been considered for the determination 
of the free parameters. The three initial free parameters 
(i.e, the two ranges and spin-orbit) were initially consid- 
ered in the fit but soon it became clear that the E/A value 
given by the polynomial fit should also be taken into ac- 
count as a free parameter in the way we discussed above. 
Out of the four, it turns out that a(E) has a very smooth 
dependency on Wls and the minimum value of cr(E) was 
always obtained for spin orbit strength values around 90 
MeV fm 5 , a value which is consistent with the value of 
130 MeV fm 5 used in Gogny D1S (the difference is related 
to the different effective masses which are 1 for BCP and 
0.7 for Gogny D1S and therefore 0.7 x 130 90). Another 
relevant observation is that the binding energy difference 
AB = B t h — -Bexp shows a linear dependence with mass 
number A (which is sometimes masked by large fluctu- 
ations at low A) with a slope that is intimately related 
to the value of E/A for the bulk. It is almost zero for 
E/A « 15.97 and is clearly different from zero and pos- 
itive for E/A = 16.03 as can be observed in Fig. O It 
turns out that the value E/A = 15.98 yields the lowest 
a(E) value. The final relevant observation is that <j(E) 
depends abruptly on the values of tql and r^jj and an 
accuracy of one part in 5 x 10 3 is required to obtain rea- 
sonable values for that quantity. It turns out that there 
are two sets of values of Tql and rou leading to reason- 
able values of cr{E), namely tql — r$u ss 0.660 fm and 
tol ~ 0.49 fm and r^u w 1.05 fm. In the latter case, 
the value of r^u is more critical than the value and 
it has to be kept at the value rou = 1.046 fm leaving 
only one free parameter to play, namely roL- Although 
the values of <j{E) at the minimum for the two choices 
rah = i"ou « 0.660 fm and rou ~ 0.49 fm and rou — 1.046 
fm are similar, the AB curve looks smoother for the first 
choice and this is the only one we will consider in the 
following. After all these considerations, and having car- 
ried out a large set of mass table calculations to validate 
our choice, we arrive to the conclusion that E/A = 15.98 
MeV, W LS = 90.5 MeV fm 5 , and r ou = r 0L = 0.659 fm 
is the best choice leading to a <r(E) value of 1.58 MeV 
for the set of 579 nuclei considered in the fit. In the 
following, as already mentioned, we will denote the new 
set of parameters as the Barcelona-Catania-Paris-Madrid 
(BCPM) functional. That the fit gives rou = rou hints to 
the fact that the ranges of the (effective) nucleon-nucleon 
forces are very similar in the different channels. Let us 
also comment on the smallness of the range parameter. 
In this respect, one should remember that even for a zero 
range force the surface energy is not zero, due to the sur- 
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face spread of the quantum single particle wave functions. 
Therefore, it only needs very little extra smearing of the 
density via a finite range force to get to the experimen- 
tal value of the surface energy. We have calculated the 
surface energy provided by the BCPM functional using 
the self-consistent Extended Thomas-Fermi approach in- 
cluding h 2 corrections (ETF-fi 2 ) [60]. According to this 
reference, the HF surface energy can be estimated by 
renormalizing the Weizsacker term in the semi-classical 
kinetic energy density. We have obtained the renormal- 
izing factor from a fit to earlier exact HF calculations 
(using BCP1 and BCP2 functionals [33]). In this way 
we estimate the HF surface energy associated with the 
BCPM functional to be 17.68 MeV. 

In panel a) of Fig. [3] we present the AB quantities 
for the 579 nuclei with measured masses as a function of 
neutron number. The first noticeable conclusion from is 
that, as expected, the agreement with experimental data 
is much better for heavy nuclei than for light ones. A 
more quantitative assessment is obtained by looking at 
a(E) for different sets of nuclei: taking into account nu- 
clei with A > 40 (536 nuclei) the a{E) for the energy 
gets reduced to 1.51 MeV and for the remaining 43 nu- 
clei with A < 40 it increases up to 2.31 MeV. When only 
nuclei with A > 80 (452) are considered a cr(E) value of 
1.34 MeV is obtained. A glance at the plot reveals that 
those figures are consistent: for heavy nuclei the fluctu- 
ations around zero are of small amplitude whereas for 
light nuclei the fluctuations lead to discrepancies as high 
as 5.5 MeV in the upper side or -5 MeV in the lower one. 
We also observe strong deviations around magic neutron 
numbers N = 82 and N = 126. In the lower panel of Fig. 
[3] we have plotted the AB for those 193 nuclei with an 
extrapolated binding energy in Audi and Wapstra compi- 
lation. As these nuclei were not considered in the fit, the 
results can be viewed as a prediction of our functional. 
We observe a very good reproduction of the binding en- 
ergies for N sa 150 but as N decreases strong deviations 
from the extrapolations are observed. For the 193 nuclei 
we obtain a(E) = 2.34 MeV which deviates quite a lot 
from the nominal value of 1.58 MeV. However, the value 
for A > 80 nuclei (146 nuclei) is 1.65 MeV which is much 
closer to our reference value. For A < 80 (47 nuclei) 
we get <t(E) — 3.76 MeV. We are tempted to conclude 
that the predictive power of BCPM for nuclei far from 
the stability line is preserved for A > 80 but is relatively 
poor for the remaining cases although part of the dis- 
crepancy could also be attributed to deficiencies in the 
extrapolation. 

In order to assess in a more quantitative way the vari- 
ance of the parameters we have performed three sets 
of calculations. In the first, the values of E/A and 
Wls have been kept and different choices for tql = tql 
have been made. The results are summarized in Table 
Mil where we observe how a change of 0.003 fm in the 
tqu = r 0L values (a 0.5 % change) modifies the cr(E) val- 
ues by more than 40 % from 1.58 to 2.16 (2.23). It is also 
worth mentioning that the a{E) value for A < 40 remains 
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Figure 3: In the upper panel, the binding energy difference 
AB = Bth — Bbx P (in MeV) is plotted with our optimal set of 
parameters (BCPM) for the 579 nuclei of Audi's AME 2003 
as a function of neutron number N. Points corresponding to 
the same isotope are connected by straight lines. In the lower 
panel, the same quantity is plotted but this time for the extra 
set of 193 nuclei in Audi's compilation with "experimental" 
values obtained from extrapolation and/or systematics. 



essentially constant as a function of tql = tqu indicating 
a high degree of randomness in AB for light nuclei that 
points to a clear deficiency of the mean field theory to de- 
scribe such light systems. On the other hand, for heavier 
nuclei with A > 40 there is a definite parabolic structure 
with the minimum centered around r$u — tql = 0.659 
fm. For A > 80 the a(E) value at the minimum gets 
reduced to 1.35 MeV in agreement with the smaller val- 
ues of AB observed for heavy nuclei. On the other hand, 
the A < 80 <j(E) values seem to favor smaller values of 
fou — ?*ol but the gain does not seem to be very sig- 
nificant. From the values in the Table we can conclude 
that the use of different values of tqu = tql for different 
regions of the mass table is marginally beneficial. 

Keeping E/A constant and equal to its optimal value 
of 15.98 MeV and varying Wls with tql = tql values 
determined as to minimize <j(E) leads to the results of 
Table IIVI There we observe that changes of the order 
of 5% in Wls lead to changes of the order of 6 % in 
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roL 


a(E) 


A > 40 


A < 40 


A > 80 


A < 80 


0.656 


2.23 


2.22 


2.35 


2.27 


2.11 


0.657 


1.91 


1.87 


2.32 


1.84 


2.11 


0.658 


1.68 


1.62 


2.31 


1.52 


2.16 


0.659 


1.58 


1.51 


2.31 


1.35 


2.23 


0.660 


1.65 


1.58 


2.33 


1.40 


2.33 


0.661 


1.85 


1.82 


2.35 


1.64 


2.45 


0.662 


2.16 


2.14 


2.40 


2.02 


2.61 





x=Wls 
j=Wls 


x=Wls 
y=roL 


x=rot/ 
y=rou 


x=roL 
y=roL 


X=ro!7 

y=r L 


dxdy 


0.195 


-166.812 


8.58xl0 4 


l.llxlO 4 


2.99 xlO 4 


a' A a(E) 
dxdij 


0.16 xlO 4 


-0.99 xlO 4 


3.72 xlO 4 


0.48 xlO 4 


1.30xl0 4 



Table III: Values of a(E) (MeV) for different values of r 0L 
(fm) in the roL = tql case. The values of E/A and Wls 
are fixed at the BCPM values of 15.98 MeV and 90.5 MeV 
fm a , respectively. Different columns correspond to different 
choices in the A values included in the evaluation of cr(E). 
For A > 40 (A > 80) 536 (452) nuclei are included. 



roL (fm) 


W LS 


a(E) 


a{E) (A > 80, 452) 


a(E) (A < 80,127) 


0.653 


85.5 


1.62 


1.38 


2.27 


0.659 


90.5 


1.58 


1.34 


2.23 


0.665 


95.5 


1.67 


1.44 


2.32 



Table IV: Values of cr(E) (in MeV) for different choices of Wls 
(in MeV fm 5 ) and at the fixed value E/A = 15.98 MeV. For 
each Wls the values roz, = roL have been fixed to minimize 
o(E). The values in parenthesis following u(E) represent the 
number of nuclei used in the evaluation of er(E). 



ct(E). This result clearly indicates that, as expected, the 
dependency of the masses on the spin-orbit strength is 
rather weak. 

Finally, the second derivatives of the a(E) with respect 
to the parameters of the fit have been evaluated. Their 
values are given in Table [V] By looking at the last row 
corresponding to the normalized variables we can read 
off how a 1% change in the corresponding value of the 
variables is going to modify the u(E) value at the min- 
imum. For Wls that 1% change will lead to a change 
of 0.16 MeV (the smallest), whereas for rou it will im- 
ply a change of 3.72 MeV (the largest). We conclude 
that the quality of our fit is more sensitive to rou than 
to tol- An eigenvalue analysis of the Hessian matrix for 
rou and rot shows that there is a nearly zero eigenvalue 
for those combinations satisfying rou — — 0.35roi- This 
implies that rou an d rou are not independent variables 
at leat around their optimal values rou = r 0L = 0.659. 
The other eigenvalue is 4.17 x 10 4 and corresponds to the 
combination rou = 2.9ro_L. 

We would like to discuss the binding energy differences 
AB from a different perspective by looking at its behav- 
ior with neutron number for each isotopic chain consid- 
ered. The quantity is plotted in Fig. 0J In each panel, 
values for the number of protons Z and the number of 
neutrons No signaling the origin of the x axis are given. 
Perpendicular lines in the AB = horizontal lines indi- 
cate the location of magic neutron numbers. From this 
plot, we conclude that for medium mass and heavy nuclei 



Table V: Values of the second derivative of a(E) (in MeV) 
with respect to the Wls (in MeV fm 5 ), rou (in fm) and roL (in 
fm) parameters. The last row is for the values corresponding 
to variables normalized to the value 1 for the parameters that 
minimize er(E). 



the AB curves are rather flat as a function of N except 
for those Z values corresponding to magic numbers (Z=50 
and 82) plus or minus two. For lighter nuclei the result 
is less certain but overall we can say that again values 
of Z of 38, 40 (not magic but a sub-shell closure) and 
Z= 28 and 30 have larger values of AB than the other Z 
values. The impact on the <j(E) can be significant. For 
instance, if nuclei with Z=48, 50, 80 and 82 are not in- 
cluded in the evaluation of a(E) its value goes down to 
1.46 for 514 nuclei. Partial <r(E) values corresponding to 
A > 80 go from 1.34 (452 nuclei) when all possible Z val- 
ues are considered to 1.09 (387 nuclei) when the Z values 
indicated above are not considered. It is also observed 
that for neutron numbers equal to magic ones plus or mi- 
nus two the binding energy difference AB always show 
peaks. The effect is reinforced when both proton and 
neutron numbers are magic, in such a way that the two 
nuclei with the largest AB values are 208 Pb (Z=82) and 
130 Cd (Z=48). Let us finally mention that the peculiari- 
ties of this Figure are also observed for other forces like 
Gogny D1S [52. 61]. The origin of these findings is uncer- 
tain but probably has to do with dynamical correlations 
associated to shape and pairing degrees of freedom. 

Another interesting piece of information is the contri- 
bution to a(E) according to the deformation of the nuclei. 
Of the 579 nuclei considered in the evaluation of <j(E) 
there are 365 with a ground state deformation parameter 
02 greater, in absolute value, than 0.1 (a value that cor- 
responds to a moderately deformed system). The <j(E) 
value for those deformed systems gets reduced to 1.08 
MeV whereas the complementary one corresponding to 
near spherical systems (214 nuclei) grows up to a value 
of 2.19 MeV. This result is in the line of recent claims 
that well deformed systems are described better by mean 
field models than the spherical ones [50(. The result is 
also consistent with the previous finding concerning the 
maximum values of AB for magic or semi-magic nuclei 
as those nuclei tend to have a nearly spherical ground 
state shape. 

To finish this subsection we would like to com- 
pare our results with those obtained with the different 
parametrizations of the Gogny force, namely D1S, DIN 
[6^ and DIM [63]. The results for a{E) and the three 
parametrizations are given in Table IVII For the three 
parametrizations the values obtained with the raw HFB 
energies as well as those including the rotational energy 



re o 
to 

B > 
tO 



re o 
95 



— 

O 
O 


c" 

r+ 

O 



N hcj 

p e 

5 to 
re 

Cfl ■ ■ 

B. > 

3 to 

B S 



B^ X 

IT 2: 
z£ 

1.1' 

£ to 
re 3 

r+ re 

B - ^ 

£ oq 

S" n & 
<" ° ft 

" 2 ™ 

re M 
.2 to 

V, 

o o 
B re 



Z,N=104, 152 

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 



Z,N=106, 154 : 

) 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 U 



'I""!""!""!""!""!""!" 



BCPM E/A=1 5.98 + FSE r ^r ,=0.659 \N LS = 90.5 




B 

o 



> 

£- 
c' 

B 
g_ 

O 
*-i 

^3 
O 
B 
g- 
n' 
B 

p" 



° % 
B - 

re 

<S. re 
B 

£» 
p p 

» 3" 

2. B 

1 a 

o' 

B 



re o 



p - 

B a 

2- Bft 

re 

> re 



B 



S < 
O p 

B B- 



§ re B 

» 8 § 

i-! re 

cn m B^ 

tr S t3 

pi ti ji 

< O p 

re p re 

P & ~ 

F B B" 

o oq re 



11 



correction (computed according to our procedure) are 
given. However, the numbers obtained in that way are 
not very fair in a comparison with BCPM's value be- 
cause in two cases (D1S and DIN) the binding energies 
of a few spherical nuclei were used in the fit, and in the 
other, the binding energy includes quadrupole zero point 
energy corrections not taken into account in our calcula- 
tions. Therefore, in order to provide with a number we 
can compare with BCPM's results we have allowed for a 
global energy shift of the binding energies with shift val- 
ues obtained as to minimize the a(E) values (also given 
in the table). Both D1S and DIN were fitted to the 
binging energies of just a few spherical nuclei and there- 
fore the impact of the rotational energy correction in the 
binding energy of deformed nuclei was not included. The 
effect of adding the rotational correction is noticeable 
bringing a(E) down from its original value of 3.48 (4.88) 
MeV to 2.15 (2.84) MeV for D1S (DIN). The additional 
phenomenological shift considered reduces the raw HFB 
a(E) value substantially in the two cases, but it barely 
affects the quantity including the rotational energy cor- 
rection in the D1S case and reduces to a surprisingly 
good value of 1.45 MeV the DIN value. On the other 
hand, DIM is fitted to essentially all the experimentally 
known binding energies and both a rotational energy cor- 
rection and a zero point energy correction associated to 
quadrupole fluctuations was considered for each nucleus 
[63. The raw HFB value given in our table is therefore 
very large but the energy shift reduces substantially the 
a(E) value down to 2 MeV. When the rotational energy 
is included a substantial reduction in a(E) is observed. 
An additional energy shift (meant in this case to mock 
up the quadrupole zero point energy correction, a quan- 
tity which is not accessible with our computing resources) 
brings the a(E) values down to 1.45 MeV. However, this 
value is still far from the 0.7 MeV given in [63]. The 
discrepancy has to be attributed to the slightly different 
rotational energy corrections and the lack of quadrupole 
zero point energies in our calculation. The numbers ob- 
tained for DIM and DIN are close to the one of BCPM 
(1.58 MeV) when the same set of nuclei and beyond 
mean field effects are considered in all the cases. This 
fact raises the expectation that a BCPM functional fit- 
ted with quadrupole zero point energies (and eventually 
taking into account the correlation energy of dynamical 
parity restoration [rJij and/or particle number symmetry) 
could lead to a much smaller value of cr(E), well below 1 
MeV. 



E. Role of the low density part in the BCPM 
functional 



One may wonder how significant is the detailed den- 
sity behavior of the microscopic input below saturation. 
For instance the low density part gives more attraction 
than most of the more phenomenological mean field ap- 
proaches. It is not the purpose of this paper to inves- 



a(E) 


D1S 


DIM 


DIN 


HFB 


3.48 


5.08 


4.88 


HFB+E fl oT 


2.15 


2.96 


2.84 


HFB + Shift 


2.53 (2.4) 


2.02 (4.7) 


2.02 (4.5) 


HFB+EflOT+Shift 


2.14 (0.2) 


1.47 (2.6) 


1.45 (2.4) 



Table VI: The a(E) values (MeV) for the three parametriza- 
tions of the Gogny force D1S [U, DIM ^ and DIN ^ and 
for different kinds of theoretical calculations. For the row 
marked as HFB the theoretical binding energy corresponds 
to the raw mean field energy, in the one denoted HFB+E^ot 
the rotational energy correction is also included. The last two 
rows correspond to the same theoretical set up as before but 
this time a global shift has been applied to the binding en- 
ergy differences as to minimize the u(E) value. The values of 
the shift parameter (in MeV) for each situation are given in 
parenthesis. 



tigate this question in detail here as it was already dis- 
cussed in an earlier paper fl8j j . However, for a first test 
of the ingredients that make up our functional, we have 
analyzed the impact of considering instead of the full 
Equation of State (EOS) a parabolic approximation to 
it. The parabolic approximation is obtained from the se- 
ries expansion of the microscopic EOS around saturation 
Ep al (p) ~ E(po) + k(p — po) 2 - It has the obvious disad- 
vantage of reaching a finite value at zero density, which 
is unphysical. In order to correct this unwanted behavior 
another quadratic approximation has also been consid- 
ered, namely i?Q ua d = + (2P 2 - It goes to zero at zero 
density, and the two parameters are adjusted to have the 
correct minimum at saturation. As a consequence of the 
limited number of parameters, it is not possible to mod- 
ify the curvature of E{p) around saturation. In Fig. [5]we 
plot the three EOS considered as a function of density. 

The finite nuclei results with the parabolic approxi- 
mation turn out to be rather poor, probably as a con- 
sequence of the behavior at zero density. On the other 
hand, the quadratic approximation results turn out to be 
quite reasonable in spite of the <j{E) value of around 10 
MeV obtained without fitting the roL — tqu range pa- 
rameters. This high value comes essentially from super- 
heavy nuclei (see below). If the range parameters are 
searched for to minimize the u{E) we obtain a much 
reduced value of 3.32 (see Fig. [5]). The binding en- 
ergy difference plot still shows the problems with the 
super-heavies mentioned above and we observe that the 
drop in the binding energy difference for neutron num- 
bers in excess of 130 is completely different from the rest 
which looks rather reasonable. The behavior of the super- 
heavies may have to do with the fact that there is an ex- 
treme balance of Coulomb, surface and shell energies and 
any slight perturbation of the full solution may lead to a 
strong failure. In conclusion, we can say that a quadratic 
fit to the EOS may work reasonably well but for optimal 
numbers, the details of the microscopic density depen- 
dence of the EOS below saturation are important. 
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Figure 5: The different equations of state considered are plot- 
ted as a function of the density p in units of the saturation 
density, po. The microscopic EOS is the one used to derive 
the BCPM functional (full line), whereas the parabolic one 
(dotted line) is the one obtained by a parabolic fit to the mi- 
croscopic EOS at saturation. Finally, the Quadratic EOS is 
given by the dashed line. In the legend the parameters of the 
quadratic EOS are shown. 



IV. PERFORMANCE OF THE NEW 
PARAMETRIZATION FOR OTHER 
QUANTITIES 

Nuclear binding energies play an important role in the 
description of atomic nuclei but there are other relevant 
observables in nuclear dynamics and therefore any func- 
tional expected to become "universal" has to perform well 
also in the description of those quantities. We will shortly 
discuss a couple of the most relevant ones. 

A. Radii 

The nuclear charge radius is a relevant observable con- 
nected to many other physical quantities. Experimen- 
tally, it can be obtained easily using electron scatter- 
ing or other complementary techniques. Theoretically, 
already at the mean field level it is possible to obtain 
reasonable values for charge radii if the proton radius of 
0.8 fm is taken into account in the final value. For this 
reason we have computed charge radii for all even-even 
nuclei considered and compared the results to the ex- 
perimental data published in Angeli's compilation (65| . 
A total of 313 even-even nuclei have a measured charge 
mean square radius. When confronted to the theoreti- 
cal predictions we get a rms of 0.0347 fm which is not 



a bad number taking into account that the charge ra- 
dius has not been considered in the fitting protocol. As a 
comparison, the same quantity for the same nuclei com- 
puted with the Gogny D1S force is 0.029 fm and with the 
more recent parametrization DIM [iH is 0.036 fm. The 
comparison is specially satisfactory with Gogny DIM as 
this interaction has been fitted using global data as with 
our functional. We notice that most of the rms value for 
BCPM comes from the heaviest nuclei where a strong de- 
viation between theory and experiment is observed. We 
also notice systematic deviations for nuclei around N=60 
and N=80, two regions of the Nuclide chart character- 
ized by the phenomenon of shape coexistence. As for 
the binding energies, the deviations strongly fluctuate in 
light nuclei with N < 40. The same pattern is also ob- 
served for the Gogny D1S results, suggesting that the 
origin for the discrepancies has to do more with a poor 
description of the ground state than with deficiencies of 
the interactions. A comparison of the two plots reveals 
that a better figure for gr could be achieved in BCPM 
if an overall displacement would be performed, i.e. the 
theoretical values systematically underestimate the ex- 
perimental values. It may be that size fluctuations (RPA 
correlations) bring the radii to their correct value. An 
increase of radii has been found in (66| . This is a task for 
the future. 



B. Quadrupole and octupole deformations 

The quadrupole deformation parameter of the ground 
state is another relevant parameter associated to low 
energy nuclear properties like 2 + excitation energies or 
B(E2) transition probabilities. The connection of the 
intrinsic deformation parameter /?2 with experimental 
B(E2) transition probability values is somehow uncer- 
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Figure 7: Differences Ar between the computed radii and the 
experimental data taken from the compilation of Angeli [65| . 
In the top panel, the results for the BCPM functional and in 
the lower one, for comparison, the results for the Gogny D1S 
force. 



tain as it relies on the strong deformation limit of an- 
gular momentum projected theories to obtain the ro- 
tational formula [67(. For this reason we have pre- 
ferred to compare BCPM's predictions with the ones of a 
well reputed, performing and predictive effective interac- 
tion. We have chosen the Gogny force, D1S, (results for 
the most recent published parametrization of the Gogny 
force, DIM, are very similar). In the upper panel of 
Fig. [8] a histogram is plotted depicting the number of 
nuclei with ground state quadrupole deformation param- 
eters fa between fain) and fa{n) + 602- The discrete 
fain) values are given by the sequence fain) = 5 fan 
with n = . . . , -2, -1, 0, 1, 2, . . . The value 8 fa = 0.0125 
has been chosen in such a way that each bin represents 
roughly 5% of the typical value of 0.25 for the quadrupole 
deformation parameter. On the lower panel we plot a 
histogram with the difference /3 2 (D1S) - /3 2 (BCPM). We 
observe how most of the 818 nuclei considered have a 
difference of less than 0.0125 in absolute value. A more 
detailed analysis of the results reveals that most of the 
discrepancies take place in the region Z s» N ps 40 which 
is widely recognized as a region of prolate-oblate shape 
coexistence and triaxiality. 

Octupole deformation is associated to the multipole 
operator of order three and breaks the parity symmetry. 
It is not as common as quadrupole deformation in the 
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Figure 8: The number of nuclei with a given quadrupole de- 
formation parameter /?2 in their ground state (in bins 0.0125 
units wide) are plotted versus 02 in the upper panel. On the 
lower one, the difference between the ground state's /?2 values 
obtained with BCPM and Gogny D1S are depicted. 



ground state of atomic nuclei as very specific combina- 
tions of orbitals around the Fermi surface have to occur 
to favor it. Known regions showing octupole deformation 
are the region around the 224 Ra (Z=88, N=136) and the 
region around 144 Ba (Z=56 and N=88). We have tested 
the octupole properties of BCPM by relaxing the parity 
symmetry in our self-consistent calculations. In this way, 
the system can end up in an octupole deformed configu- 
ration in its quest for the minimum energy. The results 
show energy gains of the order of a few hundred keV and 
up to 1 MeV in nuclei in the aforementioned regions and 
both the energy gains and values of the fa deformation 
parameters are compatible with those obtained with the 
Gogny interaction [64| ■ Overall, the octupole correlation 
energies are about a few hundred keV smaller than the 
corresponding ones obtained with Gogny D1S and DIM 
ones. Typical examples are the case of 224 Ra with oc- 
tupole correlations energies of 0.63 MeV and 1.31 MeV for 
BCPM and Gogny D1S, respectively; or the one of 144 Ba 
where the corresponding values are 0.15 MeV and 0.68 
MeV for BCPM and Gogny D1S respectively. For D1S 
there are five nuclei with mean field octupole correlation 
energies larger than 1 MeV and another nine with ener- 
gies between 1.0 and 0.5 MeV. For BCPM the numbers 
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are two with energies above 1 MeV (by just a few keV) 
and nine in the interval 1.0 to 0.5 MeV. The situation 
is similar to the case of BCP1 studied in Ref. [35j with 
the only difference that in BCP1 the mean field octupole 
correlation energy was comparable to the one of Gogny 
D1S. However, as a consequence of the larger collective 
masses, the excitation energies predicted by BCP1 were 
substantially lower than the ones of Gogny D1S and the 
experiment. It remains to be tested what would be the 
spectroscopic predictions of BCPM concerning negative 
parity excited states but the lower mean field correlation 
energies go in the right direction. 

Another testing ground for quadrupole and octupole 
properties of atomic nuclei is the shape evolution from 
the ground state to spontaneous fission. For this reason 
we have performed fission barrier calculations for a few 
selected examples and the results obtained for fission bar- 
rier heights and widths as well as mass asymmetry near 
scission (connected with fragment's mass asymmetry) are 
quite close to the results obtained with D1S. A detailed 
account of these calculations will be presented elsewhere 
and here we will focus on describing the results for a cou- 
ple of examples: one is the paradigmatic case of 240 Pu 
whose fission properties have been computed with almost 
any proposed interaction. The other is the super-heavy 
262 Sg where experimental data on spontaneous fission ex- 
ist. In both cases, the results will be compared with the 
benchmark calculations carried out with the Gogny EDF. 
For an early account of fission barrier properties with pre- 
vious versions of the BCP functional (BCP1 and BCP2) 
the reader is referred to Ref. [68j|. 

We have followed the standard procedure that consists 
in the evaluation of the HFB energy as a function of the 
mass quadrupole moment Q20 with the other multipole 
moments free to adopt any value in order to minimize the 
energy. Along this path to fission we compute the collec- 
tive inertias i?(Q2o) required to evaluate the spontaneous 
fission half life in the WKB approximation (see [68[ for 
details). In Figs. 191 and [TU1 results are presented for 240 Pu 
and 262 Sg, respectively. In the lower panels the HFB en- 
ergy as a function of Q20 is depicted for both BCPM and 
Gogny D1S cases. In the two considered nuclei the shape 
of the energy landscape looks rather similar regardless of 
the interaction considered. The maximum and minimum 
are located roughly at the same position and it is only 
the height of the barriers which changes with the interac- 
tion. Overall, Gogny D1S produces barriers higher than 
BCPM in agreement with the larger surface energy co- 
eficient of D1S. The relevance of this characteristics of 
BCPM depends upon the reduction of the height of the 
first fission barrier when triaxial shapes are included in 
the calculation. In the case of the Gogny forces, the re- 
duction can be as large as two or three MeV but we do 
not have a hint for BCPM yet. Work is in progress to 
explore the characteristics of BCPM regarding triaxial 
shapes. In panel c) of Figs. [5] and [TU1 the octupole and 
hexadecapole moments are depicted as a function of Q20 
for both interactions. The results for Gogny D1S and 



BCPM lie one on top of each other and can not be dis- 
tinguished in the plot. We observe that in 240 Pu octupole 
deformation develops in its way to fission pointing to a 
dominant mass asymmetric fission mode. On the other 
hand, 262 Sg remains reflection symmetric along the whole 
fission path and therefore symmetric fission is expected 
to be the dominant fission mode in this nucleus. In panel 
b) of both figures the particle-particle pairing energies are 
given both for protons (dashed) and neutrons (full) for 
the two interactions. Again, the agreement in the evolu- 
tion of these quantities with Q20 for the two interactions 
considered is very good but BCPM produces smaller pp 
pairing energies suggesting weaker pairing correlations. 
The collective inertias are quantities very sensitive to the 
amount of pairing correlations as they strongly depend 
on the excitation energies of the low lying two quasi- 
particle excitations. In panel d) of Figs. [9] and [10] it 
is observed that the collective inertias of BCPM are be- 
tween 2-3 times larger than the ones of Gogny D1S. The 
impact of the larger inertias (see [68] for a thorough dis- 
cussion of this issue) on the spontaneous fission half lives 
is to increase their values. This effect is the opposite of 
having lower barrier heights and the summed impact de- 
pends on the specific values of the quantities entering the 
"collective action". If the spontaneous fission half lives 
are computed in the WKB approximation we obtain for 
the 240 Pu case the values *i/ 2 (sf) = 6.2 x 10 18 s for D1S 
and ti/ 2 (sf) = 1.2 x 10 28 s for BCPM, to be compared 
to the experimental value of 3.5 x 10 18 s. For the 262 Sg 
case we obtain the values ti/2(sf) = 2.96s for D1S and 
ti/2 (sf) = 4.2s for BCPM, to be compared to the exper- 
imental value of 6.9ms. From the above results we con- 
clude that BCPM tend to yield larger values of t 1 / 2 (si) 
than Gogny D1S and those values can be several orders of 
magnitude larger than Gogny D1S. This is a direct con- 
sequence of the larger collective inertias and therefore a 
direct consequence of the reduced pairing correlations. 
This results could lead us to the conclusion that BCPM 
should include stronger pairing correlations if we want 
to improve on spontaneous half lives of fission. However, 
we have to be aware of the impact on small changes in 
the parameters entering the WKB formula. For instance, 
the ground state energy usually incorporates some sort of 
zero point energy correction (denoted eo hi the literature) 
that strongly influences the WKB values of t 1 / 2 (si). Its 
value has to be connected to the quantum effects associ- 
ated to quadrupole motion but its exact meaning is still 
uncertain. As an example of the sensitivity of t 1 / 2 (si) to 
this parameter, let us mention that if we use a value of 
eo=2.5 MeV instead of the value eo=1.5 MeV used origi- 
nally we will obtain half lives which are 6 to 12 orders of 
magnitude smaller. In systematic fission calculations [69] 
the value of eo is fine tuned for each isotopic chain consid- 
ered and therefore absolute values for some choice of eo as 
given here should not be considered seriously for a com- 
parison with experiment. Finally, the agreement between 
BCPM and Gogny D1S in the values of Q30 along the fis- 
sion path indicates that both interactions will produce, 
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Figure 9: (Color online) Fission properties of 240 Pu obtained 
with the Gogny D1S interaction (blue) and the BCPM func- 
tional (black). In the bottom panel the HFB energy is de- 
picted as a function of the mass quadrupole moment from the 
spherical configuration up to a very elongated configuration. 
In the other three panels, quantities relevant for the fission 
half life are given: from bottom to top, the particle-particle 
correlation energies for protons (dashed lines) and neutrons 
(full lines); the octupole and hexadecapole moments; and in 
the upper panel, the collective mass for quadrupole motion. 



at least in a static framework, the same fission fragment 
mass distributions. 



V. EXCITATION PROPERTIES OF THE BCP 
ENERGY DENSITY FUNCTIONAL 

Although the BCPM energy density functional is basi- 
cally build up to deal with some properties of the ground 
state such as energies and radii, it is however interesting 
to check if this functional can also provide useful infor- 
mation about some excited nuclear states, as for example 
the isoscalar giant monopole and quadrupole resonances, 
ISGMR and ISGQR, respectively. Giant resonances can 
be understood as small amplitude oscillations of nuclei 
that are the response to an external field generated by 
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Figure 10: Fission properties of the super-heavy nucleus 
262 Sg. See Fig. |9]for the figure caption. 



hadronic or electromagnetic probes. A useful theoretical 
framework for describing these oscillations is the RPA 
[BlT | that allows to obtain the particle-hole strength func- 
tion S(E) which describes the response of the nuclei. 
However, if the strength is concentrated in a narrow re- 
gion of the energy spectrum, as it is for example the case 
of the monopole and quadrupole oscillations in medium 
and heavy stable nuclei, the so called sum rule approach 
(4l| is a useful tool that allows to estimate the energy 
of the giant resonances without performing the full RPA 
calculations. The sum rules are defined as the moments 
of the strength function as 



E k S{E)E 



(12) 



The calculation of the sum rules is only easy to handle 
in few particular cases. Examples are just the ISGMR 
and ISGQR ones. In these cases the sum rules needed 
to estimate the average excitation energy of these reso- 
nances can be simplified. Of particular interest are the 
cubic-energy, the energy weighted and the inverse energy 
weighted sum rules, 1713, mi and m_i, respectively. Once 
they are determined one can obtain two estimates of the 
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Table VII: Theoretical E3 and E\ estimates of the average 
excitation energy of the GMR including pairing correlations. 
The E3 estimate of the GQR, also including pairing, is also 
displayed. The experimental energy of the centroid and the 
corresponding error for the GMR and GQR are also given. 



Nucleus 


E 3 (M) 


Ei{M) 




Exp(M) 


Exp(Q) 


90 Zr 


19.06 


18.32 


13.34 


17.81 


± 





32 


14.30 


± 





40 


144 Sm 


16.44 


15.62 


11.45 


15.40 


± 





30 


12.78 


± 





30 


208p b 


14.49 


13.84 


10.16 


13.96 


± 





20 


10.89 


± 





30 


112 Sn 


17.75 


16.83 


12.36 


16.1 


± 





1 


13.4 


± 





1 


114 Sn 


17.64 


16.75 


12.28 


15.9 


± 





1 


13.2 


± 





1 


116 Sn 


17.53 


16.66 


12.21 


15.8 


± 





1 


13.1 


± 





1 


118 Sn 


17.41 


16.55 


12.15 


15.6 


± 





1 


13.1 


± 





1 


120 Sn 


17.29 


16.43 


12.09 


15.4 


± 





2 


12.9 


± 





1 


122 Sn 


17.18 


16.32 


12.04 


15.0 


± 





2 


12.8 


± 





1 


124 Sn 


17.06 


16.21 


12.44 


14.8 


± 





2 


12.6 


± 





1 


106 Cd 


18.09 


17.07 


12.70 


16.50 


± 





19 










110 Cd 


17.85 


16.97 


12.49 


16.09 


± 





15 


13.13 


± 





66 


112 Cd 


17.74 


16.83 


12.38 


15.72 


± 





10 










114 Cd 


17.59 


16.73 


12.29 


15.59 


± 





20 










116 Cd 


17.44 


16.52 


12.19 


15.40 


± 





12 


12.50 


± 





66 



excitation energy of the ISGMR and ISGQR as 

£ 3= /™I and £ 1 = ./^r. (13) 
y mi V m_i 

The basic theory on the sum rule approach and details 
about the scaling method and constrained HF calcula- 
tions, which allow to determine the 7713 and m_i sum 
rules, respectively, is given in the Appendix. In the fol- 
lowing we will refer to the average energies provided by 
the scaling method £3 ([35]) and constrained HF calcula- 
tions E\ (|3"§|) as scaled and constrained energies, respec- 
tively. 

A. Numerical results 

Experimental information about the excitation energy 
of isoscalar giant resonances in medium and heavy nuclei 
is obtained through inelastic scattering of a particles of 
few hundreds of MeV measured at extremely forward an- 
gles. The multipole decomposition analysis allows to ex- 
tract from the scattering data the strength associated to 
different multipolarities, in particular those correspond- 
ing to monopole and quadrupole. From these strengths, 
the position and width of the ISGMR and ISGQR peaks 
can be determined [70l - [73j . 

We analyze first the excitation energies of ISGMR 
(Eqmr) of several nuclei for which the experimental val- 
ues are known [70l - [73l |. We compare these values with 
the theoretical predictions E3 and E\ provided by the 



Table VIII: E3 and E± estimates of the average excitation 
energy of the GMR without including pairing correlations. 



Nucleus 


E 3 (M) 


Ei(M) 


QO rv 

au Zr 


19.07 


18.36 


144 

Sm 


16.46 


15.77 


us Pb 


14.49 


13.84 


Sn 


17.81 


17.07 


114o 

Sn 


17.64 


17.71 


1160 
°Sn 


17.56 


16.66 


118a 

Sn 


17.42 


16.55 


120c 

Sn 


17.30 


16.44 


122 Sn 


17.19 


16.33 


124 Sn 


17.08 


16.21 


106 Cd 


18.10 


17.36 


110 Cd 


17.86 


17.11 


112 Cd 


17.74 


16.98 


114 Cd 


17.59 


16.73 


116 Cd 


17.44 


16.52 



BCPM energy density functional. The relevant informa- 
tion about these excitations energies is collected in Tables 
I VIII and I Villi As far as we are including in our study tin 
and cadmium (neglecting the effects of the small deforma- 
tion) isotopes and 144 Sm, for which pairing correlations 
are essential to describe accurately the ground state, we 
discuss first the role of pairing in the sum rule approach 
used to estimate Eqmr- 

The generalization of RPA to the case of nuclei with 
open shells when pairing correlations are active, is the 
quasi-particle RPA (QRPA). The so-called dielectric the- 
orem, which justifies constrained calculations within 
QRPA to obtain the m_i sum rule, has been proven re- 
cently by Colo and collaborators in (74j. Also Khan et al 
have shown the validity of the Thouless theorem for the 
energy weighted sum rule in the case of self-consistent 
QRPA based on HFB [7|| . As far as a simple generaliza- 
tion of the Thouless theorem allows to obtain the 7713 sum 
rule [4l| , it can be expected that this also will be the case 
in the QRPA framework. In addition we have also numer- 
ically checked that the virial theorem, i.e. the stability 
against scale transformations of the wave- functions p4|) , 
is also fulfilled in the HFB approximation. From this dis- 
cussion we may conclude that the sum rule approach can 
also be applied in the QRPA case. 

In Table IVHl we display the theoretical scaled and con- 
strained estimates (including pairing) of Eqmr. of the 
nuclei 144 Sm, 208 Pb, 90 Zr, 112 - 124 Sn, and i06,iio-iie Cd 
together with the corresponding experimental values (70| - 
|73| . These estimates are given by £3 and E\, respec- 
tively, defined previously, which, as mentioned, fulfill 
E3 > Ei (see Appendix). The escape width (l4Tfl) lies 
in a window of 2-3 MeV, pointing out that the monopole 
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Figure 11: Excitation energies of the GMR of 208 Pb and 90 Zr 
estimated with the scaling method and HF constrained cal- 
culations for different mean field models as a function of the 
parameter M defined in the text (upper panel). Excitation 
energies of the GMR of 208 Pb and 90 Zr estimated with the 
scaling approach as a function of the square root of the in- 
compressibility modulus (lower panel). 



strength of the considered nuclei is mainly concentrated 
around well defined peaks. The influence of pairing cor- 
relations on Eg mr has been discussed in earlier liter- 
ature (see [1^, uM an d references therein) . We have re- 
peated the scaled and constrained estimates of the Eg mr 
switching off pairing and the corresponding results are re- 
ported in Table Eini Comparing Tables EI and EIID we 
notice that, when pairing correlations are taken into ac- 
count (Table fVlT|) . the E\ and E$ values in open shell 
nuclei slightly decrease as compared with the values ob- 
tained without pairing (Table fVIII)) which in agreement 
with the results of Refs. [3^, l76l). The effect of pairing 



Figure 12: Density dependent incompressibility calculated us- 
ing different mean field models as a function of the density in 
the ranges 0.5 < p/po < 1-25 (upper panel) and 0.6 < p/po < 
0.8 (lower panel). 



correlations is, in general, more important in the con- 
strained calculations than in the scaled ones, especially 
in half-filled major shell nuclei as it is the case of 112 Sn, 
114 Sn, 116 Sn, 106 Cd, 110 Cd, 112 Cd. When the occupancy 
of the shell increases, as it happens in 118 Sn, 120 Sn, 122 Sn, 
124 Sn, 114 Cd and 116 Cd, the influence of pairing is actu- 
ally very small as it can be seen in Table I Villi 

In order to compare with the experiment, we will use 
as representative the constrained energy Ei, which esti- 
mates the energy of the centroid in a very precise way 
at least in the case of the nucleus in 208 Pb [nl [zl] (see 
also Table IV of [zi| in this respect). From Table |7nl 
we see that the experimental Egmr m 208 Pb is very well 
reproduced by Ei obtained using the BCPM energy den- 
sity functional. The experimental excitation energies of 
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90 Zr and 144 Sm are slightly overestimated by the theo- 
retical constrained calculation. However, our theoretical 
calculation are less successful in the case of the cadmium 
and tin isotopes. The predicted E\ values are shifted 
by around 1 MeV up as compared with the experimen- 
tal results [72L [73l|. This is a general tendency shown by 
mean- field models of different nature (79l . |8Q |. Several 
attempts to explain this behavior have been proposed in 
the past such as the density dependence of the incom- 
pressibility in neutron rich matter [8l| or different types 
of pairing interactions (3^ | . However, no clear explana- 
tion has been found to date and the question why Sn 
(and Cd) isotopes are so soft, and, therefore, it is still an 
open problem. It may be that nonlinear effects coming 
from particle vibration coupling play a role in such nuclei 
as Sn isotopes, since we reproduce the stiff doubly magic 
nucleus of 208 Pb and not the semi- magic ones as, e.g., 
the Sn isotopes. 

From the previous discussion it seems that the Eqmr 
estimated with the BCPM functional are somewhat high 
as compared with the experiment in spite of its relatively 
low incompressibility modulus Kq=21Q MeV (see Table 
HI)) . The fact that effective mean field models with dif- 
ferent values of Kq can give similar values of Eqmr hi 
208 Pb has actually been a general puzzle in the past (see 
pT7l . [79l | and references therein). The reason is the propor- 
tionality between Eqmr and the square root of Kq found 
by Blaizot [82] for Gogny forces. For example RMF ap- 
proaches systematically lead to Eqmr comparable to the 
ones obtained with non-relativistic theories like Skyrme 
or Gogny, in spite of the fact that Kq values are appre- 
ciably larger in relativistic than in non relativistic ap- 
proaches. Some time ago Colo et al. [z3| showed that 
one can construct Skyrme functionals that give reason- 
able values for Eqmr but with relatively high values of 
Kq i.e. values which approach the ones of RMF theories. 
As discussed in [77| , the underlying reason of this fact can 
be attributed to the different density dependence of the 
symmetry energy at, and around, saturation. In the same 
reference ffl\ it is concluded that within non-relativistic 
models there is not a unique relation between the value 
of Kq associated to a given model and the energy and 
the Egmr predicted by the same model. 

To clarify this latest problem a very recent proposal of 
Khan, Margueron and Vidaha [HJ may be useful. These 
authors analyzed a relatively large set of Skyrme, Gogny 
and RMF interactions finding that the density dependent 
incompressibility defined as 



dp 2 
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P(P) 



(14) 



where H is the energy density and P(p) the pressure of 
the symmetric nuclear matter, shows a crossing point at 
a density of about 75% of the saturation density. This 
crossing density corresponds roughly to the average den- 



sity of a heavy nucleus defined as 

1 



Pav 



A 



p (r)dr. 



(15) 



The authors claim that the excitation energy of finite 
nuclei is more properly correlated with the slope of in- 
compressibility K(p) computed at the crossing density 
than with the incompressibility modulus computed at 
saturation density Kq. This scenario is similar to the 
one found for the nuclear symmetry energy in bulk mat- 
ter, which also shows a crossing point for a similar value 
of the density (~0.11 fm -3 ) [84j. In this case it is also 
found [H], [85| that the slope of the symmetry energy at 
the crossing density is linearly correlated with the neu- 
tron skin of a heavy nucleus, e.g. 208 Pb. In general, the 
existence of a crossing point is actually related to the fact 
that the parameters of successful mean- field models have 
been fitted to reproduce properties of finite nuclei whose 
densities are, on the average, smaller than the saturation 
density. Therefore, if an observable measured in finite 
nuclei is related to some property of nuclear matter, the 
correlation is better described at the average (crossing) 
density than at the saturation one 18311 . 

Following the suggestion of Ref. [83], we have analyzed 
the behavior of Egmr of 208 Pb and 90 Zr calculated using 
different mean field models as a function of the parameter 
M defined as 



M 



, dK(p) t 



(16) 



which is directly related to the slope of the density de- 
pendent incompressibility K(p) at the crossing density 
Pav In our study in addition of the BCPM energy den- 
sity functional, we have also considered the Dl and D1S 
Gogny forces, the SkM* and T6 Skyrme interactions and 
the RMF parametrizations NL1, NL3 and NLSH. Using 
these models the correlation between the scaled and con- 
strained estimates of Eqmr's and the parameter M are 
displayed in Figure Qj] for the two aforementioned nu- 
clei. The data of the Eqmr predicted by Gogny forces 
and RMF parametrizations are taken from Refs. [86] and 
[87| respectively. 

From this figure it can be seen that the Egmr fol- 
low a linear trend with the parameter M (the correlation 
parameter is r=0.985 in all the considered cases) with 
high accuracy, in agreement with the suggestion of (83| . 
Consequently, and as it is claimed in (83|, the density 
dependence of K(p) can be constrained by the Eqmr 
measured in finite nuclei. The same correlation is also 
fulfilled by the nucleus 90 Zr with M evaluated at the 
same average density. This fact points out that the de- 
pendence on the mass number of the average density is 
actually rather weak for medium and heavy nuclei |83j | . 
From the lower panel of this figure we also see explicitly 
that considering globally mean field models of different 
nature, the linearity between Egmr and the square root 
of the incompressibility modulus at saturation [82] is lost 
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in some cases, in agreement with the conclusion of Ref. 

To investigate this aspect in more detail, we plot in 
Figure [12] the bulk incompressibility (fl4|) as a function 
of the density in the upper panel. In the lower panel we 
show, magnified, the relevant region of the upper panel. 
The density dependent incompressibility K{p) computed 
with Skyrme and Gogny forces and RMF takes for den- 
sities close to 0.72p/p similar values around 40 MeV. 
However, the BCPM functional predict relatively larger 
values of K(p) at this average density. The reason of 
that lies, probably, in the fact that K(p) grows more 
linearly with the density than the other considered mod- 
els. As far as the slope of K(p) has a value consistent 
with EqmRi as it can be seen in Figure [IT] the value of 
K(p) at 0.72 p/po is larger than the value predicted by 
the other considered mean field models. From the anal- 
ysis of Figure [T^] one derives two important conclusions. 
First, the crossing density actually disappears when one 
considers mean field models of different nature. What 
is actually relevant, is to compute the slope of K (p) at 
the average density. The same happens in the analysis of 
the symmetry energy using different mean field models 
[85| . Second, the small differences of K(p) at the aver- 
age density may explain why \/Kq and Eqmr are not 
always linearly correlated. This happens, e.g., with the 
NL1 parametrization with 1^00=211 MeV, a value simi- 
lar to the one of SkM* force but predicting Eqmr clearly 
smaller, or the case of BCPM discussed previously. 

In what concerns the ISGQR, we are well aware of 
the fact that a theory with m* = m, as is the case 
of the BCPM functional, underestimates the collective 
quadrupole excitation energy. In Table IVIII we display 
the scaled energies of the quadrupole excitation com- 
puted with BCPM for several nuclei together with the 
available experimental values [7014721 |88| . Insofar as the 
focus with a KS-DFT approach is on ground state prop- 
erties, we may not worry about the failure in the estimate 
of the ISGQR. In general, for the description of excita- 
tion energies, the KS-DFT should be generalized to in- 
clude non-localities. Attempts in this direction exist in 
condensed matter. It may be worth to look into that in 
the future also in the nuclear context. 



VI. CONCLUSIONS AND DISCUSSION 

In this work we further elaborated on the previously 
established Kohn-Sham DFT functional, BCP (l|, based 
on advanced microscopic nuclear and neutron matter cal- 
culations. The BCP contained by large 5 adjustable 
parameters, three for finite size, one for spin-orbit and 
one for bulk. We reduced this number by approximately 
a factor of two coming up with a KS-DFT functional 
containing only three adjustable parameters without de- 
teriorating the excellent results. One of these parame- 
ters, the strength of the spin-orbit force, has the usual 
value, the results being quite insensitive to its variation 



within a relatively large margin. We would like to call 
this a rather weak, i.e. almost predetermined param- 
eter (in addition it may be possible to extract it from 
the G-matrix [Isl). The other two parameters, on the 
contrary, are fixed within extremely tight margins of the 
order of 10 -4 . They concern the bulk energy per par- 
ticle, E/A, on the one hand and the surface energy on 
the other hand. Indeed in the polynomial fit to the mi- 
croscopically determined EOS, a single fine tuning pa- 
rameter had to be introduced in order to get a minimal 
rms value for nuclear masses as well ("IS R drift free (i.e. 
scatter around zero) difference of theoretical and exper- 
imental binding energies as a function of mass number. 
This parameter turned out to pin down the binding en- 
ergy to the value, E/A = 15.98 MeV. We emphasize that 
the fourth digit is significant here. For the second pa- 
rameter, we made the simplifying ansatz to introduce 
a finite range Gaussian convoluting the two densities in 
the p 2 terms of nuclear and neutron matter fits, that is 
p 2 — s- J d 3 rid 3 r 2 p(ri)e( ri_r2 ' / r ° p(v 2 )- The pre-factor is 
determined from the bulk, so that the only free parame- 
ter is the range ro which evidently is responsible for the 
surface energy. The minimization of the rms of masses 
(see main text for details) again is determined within the 
very narrow margin of 10~ 4 and turned out to have the 
value ro = 0.659 fm. With this, a rms value for the 
masses of 1.58 MeV for 579 nuclei and of 0.0347 fm for 
the radii with 313 nuclei is obtained. Those values are as 
good as, e.g., the ones of the Gogny D1S and DIM forces 
for the same set of nuclei. This very encouraging result 
calls for some comments. It is, indeed, quite surprising 
but on the other hand very satisfying, that the two basic 
nuclear quantities, energy per particle of the bulk and 
surface energy essentially suffice to determine a KS-DFT 
functional with excellent ground state properties of nu- 
clei. We should point out that the KS-DFT approach 
is tailored for the ground state and that a time depen- 
dent extension is not evident. We nevertheless tried to 
determine at least the GMR energies with the sum rule 
technique and found that the results are in the same ball 
park as the ones of other Skyrme, relativistic, or Gogny 
type approaches. On the other hand the energies of the 
GQR seem to be underestimated systematically by about 
1 MeV. This is not so surprising, since a particularity of 
the KS-DFT approach is that the bare mass is used while 
it is known that the GQR is sensitive to the value of the 
effective mass. 

Further developments concern eventually the three 
body force which may become better under control in 
the future. A first step in this direction could be to shift 
the one parameter which serves to fine tune the present 
microscopic EOS to a more deeper level, i.e. to an addi- 
tional parameter in the three body force. This eventually 
is a more satisfying knob to turn than the adjustment of 
the EOS at the end of the calculation. Whether such 
a procedure is possible remains to be seen. A further 
point which needs to be improved is pairing. At the 
moment we treat it at a rather elementary level with 
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a zero range force and cut off. A finite range force of 
the Gogny type in the pairing channel would certainly 
be preferable. Promising developments in this direction 
have recently been published [90]. At the end let us men- 
tion that our theory fulfills all the eleven criteria put for- 
ward in besides one. Our functional does not agree 
with the constraint concerning the value of the volume 
isospin incompressibility coefficient that appears in the 
leptodermous expansion of the finite nucleus incompress- 
ibility [i^. This quantity is defined as 

K T ,v = {Ksym ~ 6L + ^-L) . (17) 

The BCPM functional predicts K T ,^-195.2 MeV which 
is outside of the range proposed in [43||, -760 < K TV < - 
372 MeV. However, two comments are in order. First the 
extraction of this coefficient from the experimental exci- 
tation energies of the GMR in finite nuclei (43l . l72l [73j is 
not a simple task at all and may give erroneous estimates 
of the coefficients of the leptodermous expansion of the 
finite nucleus incompressibility as it has been discussed in 
the past (see for example Refs. dHHH). Second, there 
are other estimates of the K T _ V coefficients available in 
the literature, as for instance in Ref. (92|, that predict 
-370 ± 120 MeV which embraces the prediction of the 



BCPM functional. In general, the criteria and their lim- 
itations in [43] can certainly be subjected to debate. For 
example there is no reason to exclude functionals with 
to* = to as we have demonstrated in this work. The fact 
that our results comply with the criteria of [43j is an in- 
dication that nuclear density functionals should be based 
as much as possible on a solid microscopic input. In this 
respect it is to be noted that out of the five functionals 
which survived the criteria in [43j], three are also based 
on a microscopic input. 
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VII. APPENDIX: SUM RULE APPROACH TO 
GIANT MONOPOLE AND QUADRUPOLE 
RESONANCES 

In this Appendix we briefly summarize the basic the- 
oretical aspects of the sum rule approach to obtain mo- 
ments of the RPA strength functions. For more details 
we refer the reader to Refs. 

In the sum rule approach it is enough to know few 
low energy weighted moments of S(E) (the so-called sum 
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rules) to estimate the excitation energies of the giant res- 
onances. The sum rules are defined as 



J2 E n\(n\Q\0)\ 2 , 



(18) 



where Q is the operator associated to the excitation, |0 > 
and 1 7i > are the exact ground and excited states respec- 
tively and E n the excitation energies. If k is an odd 
integer, the sum rules mu can be computed as the expec- 
tation values in the ground state |0) of some commutators 
which involve the hamiltonian H and the operator (J (as- 
sumed to be a hermitian and one-body operator) [4J| , as 
for example 



and 



ii = {0\[Q,[H,Q]]\0) (19) 



m a = {0\[[Q,H],[H,[H,Q]]]\0), (20) 



from where an average excitation energy can be estimated 
as 



(21) 



The full calculation of the sum rules is still a compli- 
cated task because the exact ground-state wavefunction 
is, in general, unknown. However, if the moments are 
computed within the lplh RPA, it is possible to replace 
the exact ground-state wavefunction by the uncorrelated 
HF one to obtain the sum rules dTHJ) and [___]. In 
spite of all these simplifications, the calculation of the 
sum rules is only easy to handle in few particular cases, 
namely the ISGMR and ISGQR ones. In these cases 
the sum rules needed to estimate the average excitation 
energy of these resonances can be simplified as we will 
discuss below. 

For the isoscalar monopole and quadrupole oscillations 
the corresponding excitation operators are taken as 



i=i 



Qq 



»=i 



(rf - 3*?). 



(22) 



The underlying effective interaction associated to the 
BCPM energy density functional commutes with the op- 
erator Q and, therefore, the only contribution to the com- 
mutator [Q, [H, Q]] in eq. (TT9")) comes from the kinetic 
energy operator. Consequently (4l|, [5l[ 



2ft' 



mi 



A(r 2 ) (L = 0) mi 



4ft 2 



A(r 2 ) (L = 2), 
in in 

(23) 

where the expectation value of the operator r 2 is calcu- 
lated with the HF wave functions. 

The direct evaluation of the commutators entering in 
Eq. (|20p to compute the 7773 sum rule is a rather cum- 
bersome task, that in the case of the monopole and 
quadrupole oscillations can be avoided with the help 



the so-called scaling method [_l| . We start from the 
scaled ground-state wave function that in the case of the 
monopole reads 



4>x =A 3 / 2 o (Ax,Ay,Az), 



(24) 



where A is an arbitrary scaling parameter. In this case 
the m.3 sum rule can be expressed by means of the second 
derivative of the scaled ground-state energy [4l|, that is, 



m 3 



_____ 

29A2 



[($ A |ff|$ 



A=0' 



(25) 



In the scaling approach the BCPM energy density func- 
tional reads 

25(A) = T(X)+E% t (X)+E™(X)+E c (X)+E s -°(X), (26) 

where the scaled kinetic, spin-orbit and Coulomb energies 
can be found in Ref. [4JJ. The scaled bulk contribution 
can be easily obtained from the scaling of the particle 
density (p\(r) = \ 3 p(\r)). The scaling of the finite range 
contribution has been described in Refs. j8(| [87} ■ It is 
easy to show that the scaling of the Hartee term reduces 
to a renormalization of the range p of the interaction, 
or in other words that the scaled relative coordinate is 
sm,\ — \r — r'\/X. In this way, the second derivative with 
respect to the scaling parameter A, needed to compute 
the 771,3 sum rule, finally reads (87l |: 



dX 2 



U=i = 2 II ,/r ' // ''." ( '">('■'.) 



_ dv 

2s— + s z 
as 



ds 2 



In the BCPM energy density functional, the form factor 
v(s) is of Gaussian type with a range p. In this case the 
second derivative of the scaled finite range contribution 
can be recast as 186] 



d 2 E^(X) 2 d 2 E*(p) 
-\\=i = /J 



dX 2 



dfi 2 



(28) 



Finally the 777,3 for the monopole oscillation computed 
with the BCPM energy density functional can be written 
as 



,T n S 1 

m 3 (L = 0) = - — 

2 V 777 



2T + 20E S 



d 2 E°°(Xp) l 
dX 2 



A=l 



d 2 ElP( P ) 
dn 2 



(29) 



In the case of the quadrupole oscillation the scal- 
ing transformation of each single-particle wavefunction 
is given by [4l| 



P\(r) = p{Xx,Xy,z 2 /X 2 ). 



(30) 



Notice that under this transformation the volume el- 
ement is conserved because of the surface nature of the 
quadrupole oscillation. 
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As explained in detail in Ref. [83], a pure Hartree term 
scales in the quadrupole 

E% (*) = \fj drdr'p(r)p(P)v( SQA ), (31) 

where sq,\ — {s x /X, s y /X, X 2 s z ). After some algebra the 
surface contribution to the 7773 sum rule in the scaling 
approach for the quadrupole oscillation can be written 
as (see Ref.[87| for details) 



dX 2 U_1 5 



dfdr' p(r)p(r') 



-is 



dv 
ds 



ds 2 



(32) 

that in the case of the BCPM functional owing to the 
gaussian character of the form factor v can also be written 
as: 



d 2 E»(X) 



dX 2 



,d 2 E H (p) o dE»{p) 



dp 2 



2p- 



dp 



(33) 



This is due to the fact of the volume conservation in the 
quadrupole scaling, the bulk interaction part of the BCP 
energy density functional does not contribute to the 7773 
sum rule in this case. The contribution from kinetic and 
spin-orbit energies can be easily obtained from the scal- 
ing transformation of the single-particle wave functions 
assuming a spherically symmetric ground-state (4ll | . The 
Coulomb contribution can also be obtained from Eq. (|3"Tj) 
using the Coulomb form factor e 2 /s [87]. Collecting all 
these contributions, the 7713 sum rule for the quadrupole 
oscillation can finally be written as 



Assuming 77 sufficiently small, so that the perturbation 
theory holds, the variation of the expectation value of Q 
and H are directly related to the m_i moment |4lj . 



|(n|Q|0)| 2 _ 


1 


\d(QY 


1 


~d 2 (H)~ 


E n 


~2 


dp 


»7=0 ^ 


drj 2 



r)=0 

(36) 

Therefore at RPA level the average energy of the giant 
resonances can also be estimated performing constrained 
HF calculations, i.e. looking for the HF solutions of the 
constrained Hamiltonian 



H(r,) = H-pQ 



(37) 



where Q is the collective monopole (quadrupole) operator 
defined previously (see Eq. (|22|). From the constrained 
HF wave- function $(77) solution of (|3T|. the RPA m_i 
moment can be written as 



1 

777_1 = 

2 



^(m\QMv)) 



(38) 



— (<S>(r,)\H\<f>(r,)) 



r/=0 



i)=0 



from where another estimate of the average energy of the 
isoscalar giant resonances is given by 



£1 = 



7771 
777_l 



(39) 



m 3 {L = 2) = 4 



1W_ 

m 



2 r 



T 



1 



-E s 



1 

To 



f ,V - 2/i 



dp 2 



dE§(p) 
dp 



(34) 



Once the 7773 sum rule has been obtained through the 
scaling method, an estimate of the excitation energy of 
the ISGMR and ISGQR can be calculated with the help 
of Eqs. (HU and(|2"5|> 



Due to the fact that the RPA moments fulfill 
\J m^/mi > mi /mo > y/mx/m—x (4l| . the average ener- 
gies i?3 and Ei values are an upper and lower bound of 
the mean energy of the resonance and their difference is 
related with the variance of the strength function (reso- 
nance width) 



E, 



(40) 




(35) 



Let us now consider a nucleus, described by a Hamil- 
tonian H, under the action of a weak one-body field r/Q. 



In the RPA formalism, a is the escape width, which in- 
cludes the Landau damping but not the spreading width 
coming from more complicated excitation mechanism. 
Therefore, the width estimated with (|40l) can be signifi- 
cantly lower than the experimentally measured value [9~ij . 



